function [ DoubleDeflation_deflator_n,RGDP_doubledef_hat_n,RGDP_doubledef_hat_fj ] = ...
    DoubleDeflation_deflatorCES_fun_approx(P_hat_mnj,...
        pi_M0,pi_M_f0,pi_l_f0,pi_l_f1,pi_l_0,pi_l_1,...
        wage_hat,X_mnj_1,X_mnj_0,X_mjf_0,X_mjf_1,VA_n_0,VA_hat_fj,VA_hat_n,rho_f,rho_nj_repmat);

%% Preamble 

% Author: IM

% Purpose: Creating the real GDP using the double deflator instead of CPI
% inflation at the firm and sector level. Then using VA to isolate the 
% corresponding GDP deflator. 
%
% Following master notes3 Section 1.4

% Date: 15/07/2019


global alphaT rhoT lambdaT etaT shockT rho sigma sigmaT lambda eta ...
    tol maxit D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 



%% French PPI

% Initialising a matrix we will use to sum over k countries for France
% This is done for the intermediate price index P^M
log_P_hat_mnj=log(P_hat_mnj(:,france)');
log_P_hat_mnjf0_repmat = kron(log_P_hat_mnj,ones(FI_J,1));
log_P_hat_f_weta = pi_M_f0.*log_P_hat_mnjf0_repmat;
log_P_M_f_hat = sum(log_P_hat_f_weta,2);
P_M_f_hat = exp(log_P_M_f_hat);

log_P_M_f_hat_wlambda = log_P_M_f_hat.*(1-pi_l_f0);
% Taking wage to power of 1-lambda and multiplying by labor share
log_w_hat_f = log(wage_hat(france));
log_w_hat_f_wlambda = pi_l_f0.*repmat(log_w_hat_f,[FI_J 1]);

% Create total cost: F*1 matrix
log_b_hat_f = (log_w_hat_f_wlambda + log_P_M_f_hat_wlambda);
b_hat_f = exp(log_b_hat_f);


% Multiplying the the expression by a technology (a) shock 

p_hat_mnjf_france = b_hat_f.*a_hat_f;

% Computing the weights entering the price index
X_mj_0=sum(X_mnj_0,2);
X_mj_0_france = X_mj_0(start(france):start(france+1)-1);

for j=1:J
    omega_mjf(start_sorted(j):start_sorted(j+1)-1) = X_mjf_0(start_sorted(j):start_sorted(j+1)-1)./X_mj_0_france(j);
end

% Recovering the PPI (equation (6))
omega_p_hat_mjf = omega_mjf' .* p_hat_mnjf_france;
for j=1:J
    ppi_france(j,1) = sum(omega_p_hat_mjf(start_sorted(j):start_sorted(j+1)-1));
end

%% Input price deflator
ppi = P_hat_mnj;
ppi(start(france):start(france+1)-1,:) = repmat(ppi_france,[1,N]);
ppi(start(france)+16:start(france+1)-1,1:france-1)=1;
ppi(start(france)+16:start(france+1)-1,france+1:J)=1;

ppi_repmat = kron(ppi,ones(1,J));
pi_M_ppi_repmat = pi_M0.*log(ppi_repmat);
log_ipi = sum(pi_M_ppi_repmat,1)';
ipi =exp(log_ipi);
ipi(ipi==0)=1;
ipi_france = ipi(start(france):start(france+1)-1); 

clear pi_M_ppi_repmat ppi_repmat omega_p_hat_mjf omega_mjf;

% Recover domestic ppi for each country
for i=1:N
    ppi_((i-1)*J+1:(i-1)*J+J,1)=ppi((i-1)*J+1:(i-1)*J+J,i);
end

%% Sectoral output in both periods
X_mj_1=sum(X_mnj_1,2);
X_mj_0=sum(X_mnj_0,2);

%% Change in real output
Q_hat_mj = X_mj_1 ./ X_mj_0 ./ ppi_;
Q_hat_mj(isnan(Q_hat_mj))=1;

%% Sectoral intermediate consumptions in both periods

% French firms
% To be consistent with the denominator of the Domar weights, value of 
% intermediate consumption needs to take profits away from Sales
I_mjf_0 = X_mjf_0 .* (1-pi_l_f0).*(rho_f-1)./rho_f;
I_mjf_1 = X_mjf_1 .* (1-pi_l_f1).*(rho_f-1)./rho_f;
I_mj_france_0 = (I_mjf_0' * firmsecD_sorted)';
I_mj_france_1 = (I_mjf_1' * firmsecD_sorted)';

% Rest of the world
I_mj_0 = sum(X_mnj_0,2).*(1-pi_l_0);
I_mj_1 = sum(X_mnj_1,2).*(1-pi_l_1);
I_mj_0(start(france):start(france+1)-1,1)=I_mj_france_0;
I_mj_1(start(france):start(france+1)-1,1)=I_mj_france_1;

% Change in intermediate consumption (real)
I_hat_mj = I_mj_1 ./ I_mj_0 ./ ipi;
I_hat_mj(isnan(I_hat_mj))=1;

% Share of intermediate consumption in output
gamma_mj_france_0 = I_mj_france_0./X_mj_0_france;
gamma_mj_0 = (1-pi_l_0).*(rho_nj_repmat-1)./rho_nj_repmat;
gamma_mj_0(gamma_mj_0>1) = 1;
gamma_mj_0(start(france):start(france+1)-1,1)=gamma_mj_france_0;


RGDP_hat_mj = Q_hat_mj-gamma_mj_0.*I_hat_mj;
%RGDP_hat_mj(isnan(RGDP_hat_mj))=1;

% Domar weights
for i=1:N
    domar_mj_0(start(i):start(i+1)-1,1) = X_mj_0(start(i):start(i+1)-1,1) ./ VA_n_0(i);
end

%% Double Deflator
domar_Q_hat_minus_gamma_I_hat = domar_mj_0.*RGDP_hat_mj;
for i=1:N
    RGDP_doubledef_hat_n(i) = sum(domar_Q_hat_minus_gamma_I_hat(start(i):start(i+1)-1));
end
RGDP_doubledef_hat_n = RGDP_doubledef_hat_n';
DoubleDeflation_deflator_n = VA_hat_n ./RGDP_doubledef_hat_n;

%% Real GDP change for firms
RGDP_doubledef_hat_fj=VA_hat_fj./DoubleDeflation_deflator_n(france);

end

